!             Program block renorm-vlowk.f90
!
!             Authors:  Morten Hjorth-Jensen
!             ADDRESS:  Dept. Physics, University Oslo, N-0316 OSLO
!             E-MAIL:   morten.hjorth-jensen@fys.uio.no
!             LANGUAGE: F90/F95  
!             LAST UPGRADE : March 2007
!
!                        
!
!                    Begin setup of similarity transformation
!
SUBROUTINE  setup_vlowk
  USE partial_waves
  USE constants
  USE relcm_gmatrix
  USE single_particle_orbits
  USE wave_functions
  USE configurations
  IMPLICIT NONE
  INTEGER :: loop, number_orbits, max_conf

  !     reserve space in memory for mesh point and h.o. wave function  arrays
  ALLOCATE ( ra (n_rel), wra (n_rel));   ALLOCATE ( krel (n_rel), wkrel (n_rel))
  ALLOCATE ( rnlr (n_rel, 0:lmax, 0:nmax) )
  ALLOCATE(coulomb_relcom(0:nmax, 0:nmax, 0:lmax))
  !     set up Coulomb interaction in relative coordinates using Brody-Jacob-Moshinsky
  coulomb_relcom = 0.0_dp  
  !  The argonne v18 interaction contains the coulomb + electromagnetic corrections 
  !  explicitely
  IF ( type_of_pot /= 'argonnev18') THEN
     IF ( coulomb_included =='coulomb') CALL coulomb_integral
  ENDIF
  !     set up quantum numbers nl and NL relative and cm system
  CALL find_spdata_relcm(number_orbits)
  relcm_sp_data%max_value=number_orbits
  CALL allocate_relcm_array(relcm_sp_data, number_orbits) 
  CALL make_configurations_relcm
  !     set up mesh points for relative coordinates
  CALL vlowk_mesh                   
  !     setup h.o. wave functions, only relative coordinates needed
  CALL nocoreho_wfunction                
  !     find all partial waves with given jmin and jmax
  CALL setup_channels          
  !
  !     Note: first we find all possible configurations in the relative system
  !     only based on h.o. quantum number nl for all partial waves
  !     and allocate array for the given partial wave
  !     all configs are set up in the call to setup_configurations_rel
  !     First we search for the number of configurations in each channel
  !
  ALLOCATE ( rel_conf%nconfs_rel(no_channels) ) 
  ALLOCATE ( rel_conf%nconfs_relmodel(no_channels) ) 
  CALL number_configurations_rel(rel_conf)
  max_conf=0
  DO loop = 1, no_channels
     IF ( max_conf <  rel_conf%nconfs_rel(loop) ) THEN
        max_conf= rel_conf%nconfs_rel(loop)
     ENDIF
  ENDDO
  ALLOCATE( rel_conf%rel_ab(no_channels, max_conf))
  ALLOCATE ( v_com_rel(no_channels,max_conf,max_conf))
  v_com_rel = 0. 
  CALL setup_configurations_rel(rel_conf)      
  ! 
  !     Then we find all possible configurations in the relative and cm system
  !     based on h.o. quantum number nlNL for all partial waves
  !     and allocate array for the given partial wave
  !     all configs are set up in the call to setup_configurations_cm
  !     First we search for the number of configurations in each channel
  !
  ALLOCATE ( relcm_conf%nconfs_relcm(no_channels) ) 
  CALL number_configurations_relcm(relcm_conf)
  max_conf=0
  DO loop = 1, no_channels
     write(6,*) 'number of configs per channel: ', loop, relcm_conf%nconfs_relcm(loop)
     IF ( max_conf <  relcm_conf%nconfs_relcm(loop) ) THEN
        max_conf= relcm_conf%nconfs_relcm(loop)
     ENDIF
  ENDDO
  ALLOCATE( relcm_conf%relcm_ab(no_channels, max_conf+max_conf))
  CALL setup_configurations_relcm(relcm_conf)      
  !
  !     The configurations defined by the relative quantum numbers nl are then
  !     used to diagonalize the deuteron in momentum space 
  !     and obtain a model space effective interaction using a similarity
  !     transformation. This is done by the function setup_com_vlowk.
  !     Note that the interaction does not depend on the CoM coordinates
  !     and thus the problem can be separated out. The final effective 
  !     interaction is stored for transformation to the lab system.
  ! 
  CALL setup_com_vlowk
  !     With the final interaction in a harmonic oscillator basis defined by
  !      < nlNLJT_ZS | Veff | n'l'N'LJT_ZS > we perform the final transformation
  !     to the lab system in the function final_vlow_labsystem
  if(outputform=='pn') then
     CALL final_vlow_labsystem
  else if (outputform=='t') then
     call final_vlow_labsystem_t
  endif
  DEALLOCATE(coulomb_relcom)
  DEALLOCATE ( ra, wra, krel, wkrel) ; DEALLOCATE ( rnlr) 
  DEALLOCATE (rel_conf%rel_ab) 
  DEALLOCATE (rel_conf%nconfs_rel ) ;   DEALLOCATE (rel_conf%nconfs_relmodel ) 
  CALL deallocate_relcm_array(relcm_sp_data)
  DEALLOCATE (relcm_conf%relcm_ab) 
  DEALLOCATE (relcm_conf%nconfs_relcm ) 

END SUBROUTINE setup_vlowk
!
!           Set up the effective interaction in the lab-frame using the
!           mtx elements in the rel-CoM frame
!
SUBROUTINE final_vlow_labsystem_t
  USE single_particle_orbits  
  USE configurations
  USE constants
  USE relcm_gmatrix
  IMPLICIT NONE
  TYPE (configuration_descriptor) :: gmatrix_configs 
  REAL(KIND = 8), ALLOCATABLE :: lab_to_relcoeff(:,:)
  INTEGER, ALLOCATABLE :: lab_to_relconf(:,:)
  INTEGER, ALLOCATABLE :: lab_to_relnumber(:)
  INTEGER ::  p_parity, ang_mom, isospin,isospin_z, max_coeff, pq_confs
  COMPLEX*16, ALLOCATABLE::  gna_p(:,:),gna_m(:,:),gna_z(:,:)
  REAL(KIND=8), ALLOCATABLE:: twobody_com(:,:), twobody_r2(:,:), twobody_p2(:,:)

  !     loop over angular momenta
  DO ang_mom=j_lab_min,j_lab_max
     !     loop over isospin
     DO isospin=0,1 
        !     loop over parity values, here positive parity is 0, negative 1
        DO p_parity=0,1           
           !     find all possible configurations, large space and model space 
           CALL  number_gmatrix_confs_t&
                (ang_mom,p_parity,isospin,gmatrix_configs)
           IF (gmatrix_configs%number_confs <= 0 ) CYCLE
           pq_confs=gmatrix_configs%number_confs
           ALLOCATE(gmatrix_configs%config_ab(pq_confs+pq_confs) )
           ALLOCATE(lab_to_relnumber(pq_confs))
           ALLOCATE(twobody_p2(pq_confs, pq_confs))
           ALLOCATE(twobody_r2(pq_confs, pq_confs))
           ALLOCATE(twobody_com(pq_confs, pq_confs))

           CALL setup_gmatrix_configurations_t (ang_mom,p_parity,isospin &
                &,gmatrix_configs)
           !     find max possible number of transformation coeffs
           if(isospin==0) then
              isospin_z=0
              CALL find_ncoeffs(isospin_z,ang_mom,gmatrix_configs,max_coeff)           
              !     allocate space for various arrays needed in transform rel-cm -> lab 
              ALLOCATE(lab_to_relconf(pq_confs,max_coeff), &
                   lab_to_relcoeff(pq_confs,max_coeff))
              !     setup transformation coefficients for oscillator basis
              !     transformations from the c.m. frame to the lab frame
              CALL mosh_transf_t(isospin,isospin_z,ang_mom,gmatrix_configs, &
                   lab_to_relcoeff, lab_to_relconf, lab_to_relnumber, &
                   max_coeff)
              ALLOCATE(gna_z(pq_confs, pq_confs))
              gna_z=0.0D0
              twobody_com = 0.D0; twobody_p2 = 0.0D0; twobody_r2 = 0.0D0
              !     Performs HO transformation from rel and cm coordinates to
              !     lab system.
              CALL vlowk_free(gna_z,twobody_com,twobody_r2, twobody_p2,lab_to_relcoeff, lab_to_relconf,&
                   lab_to_relnumber,max_coeff,gmatrix_configs,isospin_z)
              CALL vlowk_print_t(isospin,p_parity,ang_mom,gna_z,gna_p,gna_m,twobody_com,twobody_r2,twobody_p2,gmatrix_configs)
              DEALLOCATE(lab_to_relcoeff, lab_to_relconf )
              DEALLOCATE(gna_z)
           else
              ALLOCATE(gna_z(pq_confs, pq_confs))
              ALLOCATE(gna_p(pq_confs, pq_confs))
              ALLOCATE(gna_m(pq_confs, pq_confs))
              gna_z=0.0D0
              gna_p=0.0D0
              gna_m=0.0D0              
              isospin_z=1
              CALL find_ncoeffs(isospin_z,ang_mom,gmatrix_configs,max_coeff)
              ALLOCATE(lab_to_relconf(pq_confs,max_coeff), &
                   lab_to_relcoeff(pq_confs,max_coeff))
              CALL mosh_transf_t(isospin,isospin_z,ang_mom,gmatrix_configs, &
                   lab_to_relcoeff, lab_to_relconf, lab_to_relnumber, &
                   max_coeff)
              twobody_com = 0.D0; twobody_p2 = 0.0D0; twobody_r2 = 0.0D0
              CALL vlowk_free(gna_p,twobody_com,twobody_r2, twobody_p2,lab_to_relcoeff, lab_to_relconf,&
                   lab_to_relnumber,max_coeff,gmatrix_configs,isospin_z)
              DEALLOCATE(lab_to_relcoeff, lab_to_relconf )


              isospin_z=0
              CALL find_ncoeffs(isospin_z,ang_mom,gmatrix_configs,max_coeff)
              ALLOCATE(lab_to_relconf(pq_confs,max_coeff), &
                   lab_to_relcoeff(pq_confs,max_coeff))
              CALL mosh_transf_t(isospin,isospin_z,ang_mom,gmatrix_configs, &
                   lab_to_relcoeff, lab_to_relconf, lab_to_relnumber, &
                   max_coeff)
              twobody_com = 0.D0; twobody_p2 = 0.0D0; twobody_r2 = 0.0D0
              CALL vlowk_free(gna_z,twobody_com,twobody_r2, twobody_p2,lab_to_relcoeff, lab_to_relconf,&
                   lab_to_relnumber,max_coeff,gmatrix_configs,isospin_z)
              DEALLOCATE(lab_to_relcoeff, lab_to_relconf )


              isospin_z=-1
              CALL find_ncoeffs(isospin_z,ang_mom,gmatrix_configs,max_coeff)
              ALLOCATE(lab_to_relconf(pq_confs,max_coeff), &
                   lab_to_relcoeff(pq_confs,max_coeff))
              CALL mosh_transf_t(isospin,isospin_z,ang_mom,gmatrix_configs, &
                   lab_to_relcoeff, lab_to_relconf, lab_to_relnumber, &
                   max_coeff)
              twobody_com = 0.D0; twobody_p2 = 0.0D0; twobody_r2 = 0.0D0
              CALL vlowk_free(gna_m,twobody_com,twobody_r2, twobody_p2,lab_to_relcoeff, lab_to_relconf,&
                   lab_to_relnumber,max_coeff,gmatrix_configs,isospin_z)
              DEALLOCATE(lab_to_relcoeff, lab_to_relconf )


              CALL vlowk_print_t(isospin,p_parity,ang_mom,gna_z,gna_p,gna_m,twobody_com,twobody_r2, twobody_p2,gmatrix_configs)
              DEALLOCATE(gna_z)
              DEALLOCATE(gna_p)
              DEALLOCATE(gna_m)

           endif
           !      free space
           DEALLOCATE(gmatrix_configs%config_ab)
           DEALLOCATE(twobody_com)
           DEALLOCATE(twobody_r2,twobody_p2)
           DEALLOCATE (lab_to_relnumber);
        ENDDO
     ENDDO
  ENDDO

END SUBROUTINE final_vlow_labsystem_t
!
!           Set up the effective interaction in the lab-frame using the
!           mtx elements in the rel-CoM frame
!
SUBROUTINE final_vlow_labsystem
  USE single_particle_orbits  
  USE configurations
  USE constants
  USE relcm_gmatrix
  IMPLICIT NONE
  TYPE (configuration_descriptor) :: gmatrix_configs 
  REAL(KIND = 8), ALLOCATABLE :: lab_to_relcoeff(:,:)
  INTEGER, ALLOCATABLE :: lab_to_relconf(:,:)
  INTEGER, ALLOCATABLE :: lab_to_relnumber(:)
  INTEGER ::  p_parity, ang_mom, isospin_z, max_coeff, pq_confs
  COMPLEX*16, ALLOCATABLE::  gna(:,:)
  REAL(KIND=8), ALLOCATABLE:: twobody_com(:,:), twobody_r2(:,:), twobody_p2(:,:)

  !     loop over isospin projection
  DO isospin_z=itzmin,itzmax 
     !     loop over parity values, here positive parity is 0, negative 1
     DO p_parity=0,1           
        !     loop over angular momenta
        DO ang_mom=j_lab_min,j_lab_max
           !     find all possible configurations, large space and model space 
           CALL  number_gmatrix_confs&
                (ang_mom,p_parity,isospin_z,gmatrix_configs)
           IF (gmatrix_configs%number_confs <= 0 ) CYCLE
           pq_confs=gmatrix_configs%number_confs
           ALLOCATE(gmatrix_configs%config_ab(pq_confs+pq_confs) )
           CALL setup_gmatrix_configurations &
                (ang_mom,p_parity,isospin_z,gmatrix_configs)
           !     find max possible number of transformation coeffs
           CALL find_ncoeffs(isospin_z,ang_mom,gmatrix_configs,max_coeff)           
           !     allocate space for various arrays needed in transform rel-cm -> lab 
           ALLOCATE(lab_to_relconf(pq_confs,max_coeff), &
                lab_to_relcoeff(pq_confs,max_coeff))
           ALLOCATE(lab_to_relnumber(pq_confs))
           ALLOCATE(twobody_p2(pq_confs, pq_confs))
           ALLOCATE(twobody_r2(pq_confs, pq_confs))
           ALLOCATE(twobody_com(pq_confs, pq_confs))
           !     setup transformation coefficients for oscillator basis
           !     transformations from the c.m. frame to the lab frame
           CALL mosh_transf(isospin_z,ang_mom,gmatrix_configs, &
                lab_to_relcoeff, lab_to_relconf, lab_to_relnumber, &
                max_coeff)
           ALLOCATE(gna(pq_confs, pq_confs))
           gna=0.0D0
           twobody_com = 0.D0; twobody_p2 = 0.0D0; twobody_r2 = 0.0D0
           !     Performs HO transformation from rel and cm coordinates to
           !     lab system.
           CALL vlowk_free(gna,twobody_com,twobody_r2, twobody_p2,lab_to_relcoeff, lab_to_relconf,&
                lab_to_relnumber,max_coeff,gmatrix_configs,isospin_z)
           CALL vlowk_print(isospin_z,p_parity,ang_mom,gna,twobody_com,twobody_r2,twobody_p2,gmatrix_configs)
           !      free space
           DEALLOCATE(gmatrix_configs%config_ab)
           DEALLOCATE(twobody_com)
           DEALLOCATE(twobody_r2,twobody_p2)
           DEALLOCATE(lab_to_relcoeff, lab_to_relconf )
           DEALLOCATE (lab_to_relnumber); DEALLOCATE(gna)
        ENDDO
     ENDDO
  ENDDO

END SUBROUTINE final_vlow_labsystem

SUBROUTINE vlowk_print_t(it,ip,ij,gna_z,gna_p,gna_m,twobody_com,twobody_r2, twobody_p2,gmatrix_configs)
  USE constants
  USE configurations
  USE single_particle_orbits
  USE relcm_gmatrix
  USE wave_functions
  IMPLICIT NONE
  TYPE (configuration_descriptor), INTENT(IN)  :: gmatrix_configs
  INTEGER :: i,j, ijd, it, ip, ij, ia, ib, ic,id
  COMPLEX*16, DIMENSION(gmatrix_configs%number_confs, &
       gmatrix_configs%number_confs), INTENT(IN)  :: gna_z,gna_p,gna_m
  REAL(KIND=8), DIMENSION(gmatrix_configs%number_confs, &
       gmatrix_configs%number_confs), INTENT(IN) :: twobody_com,twobody_r2, twobody_p2
  REAL(KIND=8)::zero
  ijd=ij+ij
  zero=0.d0
  DO i=1,gmatrix_configs%number_confs
     ia= gmatrix_configs%config_ab(i*2-1)
     ib= gmatrix_configs%config_ab(i*2)
     ia=(ia+1)/2
     ib=(ib+1)/2
     DO j=i,gmatrix_configs%number_confs
        ic= gmatrix_configs%config_ab(j+j-1)
        id= gmatrix_configs%config_ab(j+j)
        ic=(ic+1)/2
        id=(id+1)/2
        IF (it == 0) THEN
              WRITE(7,1001) ij, it, ip, ia, ib, ic, id, REAL(gna_z(j,i)), &
                   zero,zero,twobody_com(j,i), twobody_r2(j,i), twobody_p2(j,i)       
        ELSE

              WRITE(7,1001) ij, it, ip, ia, ib, ic, id, & 
                   REAL(gna_z(j,i)), REAL(gna_p(j,i)),REAL(gna_m(j,i)),&
                   twobody_com(j,i), twobody_r2(j,i), twobody_p2(j,i)
        END IF
     ENDDO
  ENDDO
1001 FORMAT(7I4,2X,6(2X,E13.6))
END SUBROUTINE vlowk_print_t
!
!                    Print out veff-mtx
!
SUBROUTINE vlowk_print(it,ip,ij,gna,twobody_com,twobody_r2, twobody_p2,gmatrix_configs)
  USE constants
  USE configurations
  USE single_particle_orbits
  USE relcm_gmatrix
  USE wave_functions
  IMPLICIT NONE
  TYPE (configuration_descriptor), INTENT(IN)  :: gmatrix_configs
  INTEGER :: i,j, ijd, it, ip, ij, ia, ib, ic,id
  COMPLEX*16, DIMENSION(gmatrix_configs%number_confs, &
       gmatrix_configs%number_confs), INTENT(IN)  :: gna
  REAL(KIND=8), DIMENSION(gmatrix_configs%number_confs, &
       gmatrix_configs%number_confs), INTENT(IN) :: twobody_com,twobody_r2, twobody_p2
  ijd=ij+ij
  DO i=1,gmatrix_configs%number_confs
     ia= gmatrix_configs%config_ab(i*2-1)
     ib= gmatrix_configs%config_ab(i*2)
     DO j=i,gmatrix_configs%number_confs
        ic= gmatrix_configs%config_ab(j+j-1)
        id= gmatrix_configs%config_ab(j+j)
        WRITE(7,1001) it, ip, ijd, ia, ib, ic, id, REAL(gna(j,i)), &
             twobody_com(j,i), twobody_r2(j,i), twobody_p2(j,i)       
     ENDDO
  ENDDO
1001 FORMAT(7I4,5X,4(5X,E13.6))
END SUBROUTINE vlowk_print
!
!          Find veff-mtx and make transformation from CoM and
!          Rel system to lab system
!
SUBROUTINE vlowk_free(gfree,twobody_com,twobody_r2, twobody_p2,lab_to_relcoeff, lab_to_relconf,&
     lab_to_relnumber,max_coeff,gmatrix_configs,isospin_z)
  USE configurations
  USE constants
  USE single_particle_orbits
  USE partial_waves
  USE relcm_gmatrix
  IMPLICIT NONE
  TYPE (configuration_descriptor) , INTENT(IN)  :: gmatrix_configs
  INTEGER :: k2,k1,i,j,nlc1,nlc2, conf1, conf2, lcm1, lcm2, chan1, chan2, as, ls, ns,&
       max_coeff, isospin_z, ncm1, ncm2, lr1, lr2, nr1, nr2, bra, ket, state 
  COMPLEX*16, DIMENSION(gmatrix_configs%number_confs, &
       gmatrix_configs%number_confs), INTENT(INOUT) :: gfree
  REAL(KIND=8) :: dcoem, e_com, e_r2, e_rcm, e_p2, e_p2com
  REAL(KIND=8), DIMENSION(gmatrix_configs%number_confs,max_coeff), &
       INTENT(IN)  :: lab_to_relcoeff
  INTEGER, DIMENSION(gmatrix_configs%number_confs,max_coeff), &
       INTENT(IN)  :: lab_to_relconf
  INTEGER, DIMENSION(gmatrix_configs%number_confs), &
       INTENT(IN)  :: lab_to_relnumber
  REAL(KIND=8), DIMENSION(gmatrix_configs%number_confs, &
       gmatrix_configs%number_confs), INTENT(INOUT) :: twobody_com, twobody_r2, twobody_p2

  gfree = 0.
  DO k2=1,gmatrix_configs%number_confs
     DO k1=1,gmatrix_configs%number_confs
        DO j=1,lab_to_relnumber(k2)
           nlc2=lab_to_relconf(k2,j)
           conf2=nlc2/1000;  chan2=nlc2-conf2*1000
           lcm2=relcm_sp_data%lrel(relcm_conf%relcm_ab(chan2,conf2+conf2))
           ncm2=relcm_sp_data%nrel(relcm_conf%relcm_ab(chan2,conf2+conf2))
           DO i=1,lab_to_relnumber(k1)
              nlc1=lab_to_relconf(k1,i)
              conf1=nlc1/1000;  chan1=nlc1-conf1*1000
              IF ( chan1 /= chan2 ) CYCLE
              lcm1=relcm_sp_data%lrel(relcm_conf%relcm_ab(chan1,conf1+conf1))
              IF ( lcm1 /= lcm2 ) CYCLE   ! cm orbital mom must be equal
              ncm1=relcm_sp_data%nrel(relcm_conf%relcm_ab(chan1,conf1+conf1))
              lr1=relcm_sp_data%lrel(relcm_conf%relcm_ab(chan1,conf1+conf1-1))
              nr1=relcm_sp_data%nrel(relcm_conf%relcm_ab(chan1,conf1+conf1-1))
              lr2=relcm_sp_data%lrel(relcm_conf%relcm_ab(chan2,conf2+conf2-1))
              nr2=relcm_sp_data%nrel(relcm_conf%relcm_ab(chan2,conf2+conf2-1))
              ! find bra and ket  states for nl and n'l'
              e_com = 0.0D0; e_r2 = 0.0D0; e_rcm = 0.0D0; e_p2=0.0D0; e_p2com= 0.0D0
              !             explicit computation of CoM and rel part
              IF ( (ncm1== ncm2) .AND.( lr1 == lr2 )) THEN
                 IF ( nr1 == nr2-1) e_r2 = -SQRT(nr2*(nr2+lr2+0.5))
                 IF ( nr1 == nr2-1) e_p2 = SQRT(nr2*(nr2+lr2+0.5))
                 IF ( nr1 == nr2) e_r2 =  2*nr2+lr2+1.5
                 IF ( nr1 == nr2) e_p2 =  2*nr2+lr2+1.5
                 IF ( nr1 == nr2+1) e_r2 = -SQRT((nr2+1.)*(nr2+lr2+1.5))
                 IF ( nr1 == nr2+1) e_p2 = SQRT((nr2+1.)*(nr2+lr2+1.5))
              ENDIF
              IF ( (nr1 == nr2) .AND. ( lr1 == lr2 ) ) THEN
                 IF ( ncm1 == ncm2-1) e_rcm = -SQRT(ncm2*(ncm2+lcm2+0.5))
                 IF ( ncm1 == ncm2-1) e_p2com = SQRT(ncm2*(ncm2+lcm2+0.5))
                 IF ( ncm1 == ncm2) e_rcm =  2*ncm2+lcm2+1.5
                 IF ( ncm1 == ncm2) e_p2com =  2*ncm2+lcm2+1.5
                 IF ( ncm1 == ncm2+1) e_rcm = -SQRT((ncm2+1.)*(ncm2+lcm2+1.5))
                 IF ( ncm1 == ncm2+1) e_p2com = SQRT((ncm2+1.)*(ncm2+lcm2+1.5))
              ENDIF
              IF ( ( lr1 == lr2 ).AND.(ncm1== ncm2).AND.( nr1 == nr2)) THEN
                 e_com =  2*ncm1+lcm1-2*nr2-lr2
              ENDIF
              bra = 0; ket = 0
              DO state =1, rel_conf%nconfs_rel(chan1)
                 as=rel_conf%rel_ab(chan1,state)
                 ns=relcm_sp_data%nrel(as)
                 ls=relcm_sp_data%lrel(as)
                 IF ( (ns == nr1) .AND. (ls == lr1)) ket = state
                 IF ( (ns == nr2) .AND. (ls == lr2)) bra = state
              ENDDO
              dcoem=lab_to_relcoeff(k1,i)*lab_to_relcoeff(k2,j)
              !  The interaction is diagonal in N
              !  this is not the case for  the G-matrix. Note also that the
              !  CoM part plus the twobody term -m\omega^2(r_i-r_j)^2/A is 
              !  independent of N, N' and L.
              IF ( ncm1 == ncm2 ) THEN  
                 twobody_com(k1,k2)=twobody_com(k1,k2)+dcoem*e_com
                 gfree(k1,k2)= gfree(k1,k2)+dcoem*v_com_rel(chan1,bra,ket) 
                 twobody_r2(k1,k2) = twobody_r2(k1,k2)+dcoem*e_r2
              ENDIF
              ! The p_ip_j correction depends however on N and N' and n and n'
              twobody_p2(k1,k2)=twobody_p2(k1,k2)+dcoem*(e_p2com-e_p2)*0.5D0  
           ENDDO
        ENDDO
     ENDDO
  ENDDO

END SUBROUTINE vlowk_free
!
!           Set up the interaction in the relative system using
!           a similarity transformation in momentum space. Quantum numbers are 
!           < nl JT_ZS | Veff | n'l'JT_ZS >
!
SUBROUTINE setup_com_vlowk
  USE single_particle_orbits
  USE configurations
  USE constants
  USE relcm_gmatrix
  USE partial_waves
  IMPLICIT NONE
  INTEGER ::  pq_confs, loop, mspace, bra, ket
  REAL(KIND = 8), ALLOCATABLE :: v2body(:,:)

  DO loop=1,no_channels   
     pq_confs = rel_conf%nconfs_rel(loop)
     WRITE(6,*) ' Total states and model states for channel: ', loop, pq_confs
     IF ( pq_confs <= 0 ) CYCLE
     ALLOCATE(v2body(pq_confs, pq_confs))
     v2body = 0.
     !     Setup the interaction
     CALL vlowk_channel(loop,v2body)
     DO ket = 1, pq_confs
        DO bra = 1, pq_confs
           v_com_rel(loop,bra,ket) = v2body(bra,ket)
        ENDDO
     ENDDO
     !  then we set up the final CoM and relative coordinate effective interaction
     !      free space
     DEALLOCATE(v2body)
  ENDDO

END SUBROUTINE setup_com_vlowk
!
!     Obtain the bare interaction in a harmonic oscillator 
!     basis plus the kinetic energy and the Coulomb part. It contains
!     also the CoM correction to the relative coordinates. The latter depends
!     on the mass number of the nucleus
! 
SUBROUTINE vlowk_channel(i,vint)
  USE wave_functions
  USE relcm_gmatrix
  USE partial_waves
  USE configurations
  USE single_particle_orbits
  USE constants
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: i  ! loop variable over NN channels
  INTEGER :: ncoup, n,np, bra, ket, k1, k2, a, c
  INTEGER :: la, lb, jang, lim1, lim2, lim3, lim4
  REAL(KIND = 8) :: e_coulomb, vsum
  COMPLEX*16, ALLOCATABLE, DIMENSION(:,:) :: heff(:,:), vzz(:,:)
  REAL(KIND = 8), ALLOCATABLE :: vkk(:,:)
  REAL(KIND = 8),  INTENT(INOUT):: &
       vint(rel_conf%nconfs_rel(i),rel_conf%nconfs_rel(i))

  ncoup=1
  IF ( orb_lrel_max(i) == jang_rel(i) ) ncoup = 1
  IF ( orb_lrel_max(i) /= jang_rel(i) ) ncoup = 2
  ALLOCATE(vkk (ncoup*n_rel,ncoup*n_rel))
  ALLOCATE(vzz (ncoup*n_rel,ncoup*n_rel))
  vzz = 0.0D0
  jang=jang_rel(i)
  !     setup the NN interaction
  vkk = 0.0D0
  CALL nocorepotential_interface(ncoup,vkk,i)
  vzz = vkk
  WRITE(6,'(12H Channel Nr:,I3,7H l_min:,I3,7H l_max:,I3,3H J:,I3,3H S:,I3,4H Tz:,I3)') &
       i, orb_lrel_min(i), orb_lrel_max(i), jang_rel(i), spin_rel(i), iso(i)
  ALLOCATE(heff (ncoup*n_k1,ncoup*n_k1))
  heff =0.0D0
  IF ((type_of_pot == 'OPEP').OR.(type_of_pot == 'Tensorinteraction')   &
       .OR.(type_of_pot == 'LSinteraction') ) THEN
     heff = vzz
  ELSEIF (type_of_renormv =='vbare') THEN
      heff = vzz
  ELSE
     ! get similarity transformed interaction
     CALL vlowk_mtx(ncoup,vzz,heff,i)
  ENDIF
  !     make now transformation to h.o. basis in the rel and cm system
  !     loop over all cm and rel coordinate configurations
  DO bra =1, rel_conf%nconfs_rel(i)
     k1=bra
     a=rel_conf%rel_ab(i,k1)
     n=relcm_sp_data%nrel(a)
     la=relcm_sp_data%lrel(a)
     IF ((n+n+la) > nlmax) CYCLE
     DO ket =1,  rel_conf%nconfs_rel(i) 
        k2=ket
        c=rel_conf%rel_ab(i,k2)
        np=relcm_sp_data%nrel(c)        
        lb=relcm_sp_data%lrel(c)
        !  No dependence of the bare interaction upon the CoM momenta
        !  The Hamiltonian is also diagonal in L and N 
        vsum = 0.
        IF ((np+np+lb) > nlmax) CYCLE
        IF ( ncoup == 1) THEN
           lim1=1; lim2=n_k1 ; lim3=1 ; lim4=n_k1
        ELSEIF ( ncoup == 2 ) THEN
           IF ( (la == lb).AND. ( jang > la) ) THEN
              lim1=1; lim2=n_k1 ; lim3=1 ; lim4=n_k1
           ELSEIF ( (la == lb).AND. ( jang < la) ) THEN
              lim1=1+n_k1; lim2=n_k1+n_k1 ; lim3=1+n_k1 ; lim4=n_k1+n_k1
           ELSEIF ( la >  lb ) THEN
              lim1=1+n_k1; lim2=n_k1+n_k1 ; lim3=1 ; lim4=n_k1
           ELSEIF ( la <  lb ) THEN
              lim1=1; lim2=n_k1 ; lim3=1+n_k1 ; lim4=n_k1+n_k1
           ENDIF
        ENDIF
        CALL vlowk_hosc(rnlr(:,lb,np), rnlr(:,la,n),REAL(heff(lim1:lim2,lim3:lim4)),bra,ket,vsum)
        !       Here we set up the Coulomb part in an oscillator basis. 
        e_coulomb = 0.
        IF ( ( la == lb ).AND.( iso(i) == -1 ))  THEN
           e_coulomb = coulomb_relcom(n,np,la) 
        ENDIF
        !  only potential energy : coulomb + V_NN
        vint(ket,bra)=vsum+e_coulomb
     ENDDO
  ENDDO
  DEALLOCATE(heff); DEALLOCATE(vkk); DEALLOCATE(vzz)

END SUBROUTINE vlowk_channel

!
!       Compute the T(G)-mtx for the actual channel 
!
SUBROUTINE vlowk_mtx(ncoup,vzz,heff,ichan)
  USE wave_functions
  USE partial_waves
  USE constants
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: ncoup, ichan
  COMPLEX*16, DIMENSION(ncoup*n_k,ncoup*n_k), INTENT(IN) :: vzz
  COMPLEX*16, DIMENSION(ncoup*n_k1,ncoup*n_k1), INTENT(INOUT) :: heff
  INTEGER :: i, j, ntot, i1, np, nq, j1,j2, tisoz 
  COMPLEX*16, ALLOCATABLE, DIMENSION(:,:) :: ham, cvec, temp
  COMPLEX*16, ALLOCATABLE, DIMENSION(:) ::  ceig
  COMPLEX*16 :: sum, input_energy

  ! dimension of vectors and matrices
  ntot = ncoup*n_k
  ALLOCATE( ham(ntot, ntot), cvec(ntot,ntot), ceig(ntot), temp(ntot,ntot)) 
  ! setup hamiltonian to be diagonalized
  ham = DCMPLX(0.0D0,0.0D0) 
  cvec = dcmplx(0.0d0, 0.0d0)
  ceig = dcmplx(0.0d0, 0.0d0)

  CALL complex_hamiltonian(ham,vzz, ntot, ichan)
  temp = ham
  CALL vlowkdiag_exact( temp, cvec, ceig, ntot )
  !
  ! construct renormalized nn-interaction via the Lee-Suzuki sim. transform
  !
  ! size of model space     : np = n_k1*ncoup
  ! size of complement space: nq = n_k2*ncoup
  np = n_k1*ncoup; nq = n_k2*ncoup
  CALL effective_int( np, nq, ntot, ncoup, cvec, ceig, heff )   
  !
  ! subtract diagonal elements to obtain Veff, 
  ! and compare exact with effective interactions: 
  !
  tisoz = iso(ichan) 
  DO i = 1, np
     i1 = i
     IF ( i > n_k1 ) i1 = i-n_k1
     heff(i,i) = heff(i,i) - ( krel(i1) * krel(i1) )/p_mass(tisoz)
     DO j = 1, np
        j1 = j
        IF ( j > n_k1 ) j1 = j-n_k1
        heff(i,j) = heff(i,j)/SQRT( wkrel(i1)*wkrel(j1) )/krel(i1)/krel(j1)
     ENDDO
  ENDDO
  DEALLOCATE(ham, temp, cvec );   DEALLOCATE(ceig)

END SUBROUTINE vlowk_mtx
!
! complex scaled hamiltonian in momentum representation 
!
SUBROUTINE complex_hamiltonian(h,vzz,ntot,ichan)
  USE wave_functions
  USE partial_waves
  USE constants
  IMPLICIT NONE
  REAL(KIND = 8) :: delta
  INTEGER, INTENT(IN) :: ntot, ichan
  COMPLEX*16, DIMENSION(ntot,ntot), INTENT(INOUT) :: h
  COMPLEX*16, DIMENSION(ntot,ntot), INTENT(IN) :: vzz
  INTEGER :: i, j, i1, i2, tisoz
  COMPLEX*16 :: wi

  tisoz = iso(ichan) 
  h = 0.
  DO i = 1, ntot 
     i1 = i
     IF ( i > n_k ) i1 = i-n_k
     h(i,i) = ( krel(i1) * krel(i1) )/p_mass(tisoz)  + & 
          (krel(i1)) * (krel(i1)) * wkrel(i1) * Vzz(i,i) 
     DO j = 1, ntot 
        i2 = j
        IF ( j > n_k ) i2 = j-n_k
        IF (i /= j ) THEN
           h(i,j) = SQRT( wkrel(i2) * wkrel(i1) ) *  krel(i1) * krel(i2) * Vzz(i,j)   
        ENDIF
     ENDDO
  ENDDO

END SUBROUTINE complex_hamiltonian
!
! calculate effective interaction for gamow shell model calcs
!
SUBROUTINE effective_int( np, nq, ntot, ncoup, cvec, ceig, heff )   
  USE constants
  USE wave_functions
  IMPLICIT NONE
  INTEGER, INTENT(in) ::  np, nq, ntot, ncoup
  INTEGER :: npp, k1, i1,i,k2,nqq,j
  INTEGER ::  model(np), orb, orbit_no(ntot)
  COMPLEX*16, DIMENSION(ntot,ntot), INTENT(in) :: cvec
  COMPLEX*16, DIMENSION(ntot,ntot) :: cvec_temp
  COMPLEX*16, DIMENSION(ntot), INTENT(in) :: ceig
  COMPLEX*16 :: cvec_pp(np,np), cvec_qp(nq,np), cvec_model(nq,np) 
  COMPLEX*16 :: ceig_p(np), ceig_model(np)
  REAL(kind = 8),  DIMENSION(ntot) ::  temp
  REAL(kind = 8) :: cvec_max, overlap(ntot)
  COMPLEX*16 :: e_a, e_b
  COMPLEX*16, DIMENSION(np,np), INTENT(inout) :: heff
  REAL*8 :: a1,a2,b1,b2

  ! calculate P-space overlap of all eigenvectors 
  ! loop over all eigen vectors
  DO orb = 1, ntot
     overlap(orb) = 0.
     DO i = 1, n_k1
        IF ( ncoup == 1 )THEN
           overlap(orb) = overlap(orb) + ABS( cvec(i,orb) )**2
        ELSEIF( ncoup == 2 ) THEN
           overlap(orb) = overlap(orb) + ABS( cvec(i,orb) )**2 + ABS( cvec(n_k + i,orb) )**2
        END IF
     END DO
     orbit_no(orb) = orb
  END DO
  ! sort all overlaps and corresponding orbit numbers 
  CALL eigenvalue_sort(overlap,orbit_no, ntot)
  cvec_pp = 0.; cvec_qp = 0.
  ! Set up eigenvalues and eigenvectors for model space and excluded space
  IF ( ncoup == 1 ) THEN
     DO i=1, np
        ! loop over all model space coefficients of exact eigenvectors |k>
        ! 
        k1 = orbit_no(i)
        DO j = 1, n_k
           IF ( j <= n_k1 ) THEN
              cvec_pp(j,i) = cvec(j,k1)
              ceig_p(i) = ceig(k1)
           ELSEIF ( j > n_k1 ) THEN
              cvec_qp(j-n_k1,i) = cvec(j,k1)
           ENDIF
        ENDDO
     ENDDO
  END IF
  IF ( ncoup == 2 ) THEN
     DO i=1, np
        ! loop over all model space coefficients of exact eigenvectors |k>
        ! 
        k1 = orbit_no(i)
        ceig_p(i) = ceig(k1)
        DO j = 1, n_k1
           cvec_pp(j,i) = cvec(j,k1)
           cvec_pp(j+n_k1,i) = cvec(j+n_k,k1)
        ENDDO
     ENDDO
     DO i=1, np
        k1 = orbit_no(i)
        DO j = 1, n_k2
           cvec_qp(j,i) = cvec(n_k1+j,k1)
           cvec_qp(j+n_k2,i) = cvec(n_k+n_k1+j,k1)
        END DO
     END DO
  END IF
  heff = 0.
  CALL lee_suzuki2( cvec_pp, cvec_qp, ceig_p, np, nq, heff )
  !  subtract diagonal terms

END SUBROUTINE effective_int
!
! eigenvalue sort
! sort cmplx vector real(a(1))<real(a(2)) < ... < real(a(n))
!
SUBROUTINE eigenvalue_sort(A,B, n)
  IMPLICIT NONE
  INTEGER :: i, j, n
  REAL(kind = 8), DIMENSION(n), INTENT(INOUT) :: A
  INTEGER, DIMENSION(n), INTENT(inout) :: B
  REAL(kind = 8) :: temp, temp1
  REAL(kind = 8), DIMENSION(n) :: temp2
  INTEGER :: orb

  DO i = 1, n
     DO j = 1, n
        IF ( ABS( A(i) )  > ABS( A(j) ) ) THEN
           temp = A(i)
           A(i) = A(j) 
           A(j) = temp
           orb = B(i) 
           B(i) = B(j)
           B(j) = orb
        END IF
     END DO
  END DO

END SUBROUTINE eigenvalue_sort
!
! Lee Suzuki similarity transformation 
!
SUBROUTINE lee_suzuki2( cvec_pp, cvec_qp, ceig, np, nq, heff )
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: np, nq
  COMPLEX*16, DIMENSION(np,np), INTENT(IN) :: cvec_pp
  COMPLEX*16, DIMENSION(nq,np), INTENT(IN)  :: cvec_qp
  COMPLEX*16, DIMENSION(np), INTENT(IN) :: ceig
  COMPLEX*16, DIMENSION(np,np) :: temp, omega2, sim1, sim2, omega_2inv, u, u_inv, & 
       cvec_pp_inv
  COMPLEX*16, DIMENSION(nq,np) :: omega
  COMPLEX*16, DIMENSION(np) :: ceig_p, omega2_eig, eigval2
  COMPLEX*16, DIMENSION(np,np) :: eigen_vec, vl, heff_rhs
  COMPLEX*16, DIMENSION(np,np), INTENT(INOUT) ::  heff
  REAL(KIND = 8), DIMENSION(2*np) :: rwork
  COMPLEX*16, DIMENSION(10000) :: work1
  COMPLEX*16 :: d, sum1, temp1(np,np), temp2(np,nq), norm, determinant
  INTEGER :: i_p,j_p, j_q, ii, jj, k, i_q , a_p, a_pp 
  INTEGER :: i, lda, ldb, ldvl, ldvr, info, lwork, ilo , ihi
  CHARACTER*1 :: jobvl, jobvr, balanc, sense
  REAL(KIND = 8), DIMENSION(np) :: scale, rconde, rcondv
  REAL(KIND = 8) :: abnrm
  INTEGER :: ipiv(np) ,j
  REAL(KIND = 8) :: a1, a2, b1, b2

  balanc = 'n';  jobvl = 'n' ;  jobvr = 'v';  sense = 'n';  lda = np
  ldvl = 1;  ldvr = np;  lwork = 10000
  eigen_vec = 0. 
  temp1 = TRANSPOSE(cvec_pp) 
  CALL zgeev( jobvl, jobvr, np, temp1, lda, omega2_eig, vl, ldvl, eigen_vec, ldvr, &
       work1, lwork, rwork, info )
  determinant = PRODUCT(omega2_eig(:))
  !  write(6,*) 'check determinant', determinant
  ! the P->Q space transformation matrix, omega 
  cvec_pp_inv = cvec_pp
  CALL cmplxmatinv(cvec_pp_inv, np, d)
  DO i_p = 1, np
     DO i_q = 1, nq
        omega(i_q,i_p) = SUM( cvec_pp_inv(:,i_p)*cvec_qp(i_q,:) )
     ENDDO
  ENDDO
  ! non-hermitian effective interaction
  ! setup 2-p effective interaction in P-space
  heff = 0.
  DO i_p = 1, np
     DO j_p = 1, np 
        heff(i_p,j_p) = SUM( cvec_pp(i_p,:)*ceig(:)*cvec_pp(j_p,:) ) 
        sum1 = 0.
        DO k = 1, np
           DO i_q = 1, nq
              sum1 = sum1 + cvec_pp(i_p,k) * ceig(k) * cvec_qp(i_q,k) * omega(i_q,j_p) 
           ENDDO
        ENDDO
        heff(i_p,j_p) = heff(i_p,j_p) + sum1
     ENDDO
  ENDDO
  ! organizing the matrix (P(1 + omega * omega)P)
  omega2 = MATMUL(TRANSPOSE(cvec_pp_inv),cvec_pp_inv)
  ! calculate sqrt and inverse sqrt of matrix (P(1 + omega * omega)P)
  CALL sqrtmat( omega2, U, U_inv ,np ) 
  heff_rhs =  MATMUL( U, MATMUL( heff,u_inv ) ) 
  ! check if heff is symmetrized:
  heff = 0.
  heff = heff_rhs
  DO i_p = 1, np
     DO j_p = 1, np
        ! make heff manifestly symmetric
        IF ( i_p /= j_p ) THEN
           IF ( ABS(heff(i_p,j_p)- heff(j_p,i_p)) < 1.E-6 ) CYCLE 
           !           WRITE(6,*) 'sym test', heff(i_p,j_p), heff(j_p,i_p)
        ENDIF
     ENDDO
  ENDDO
  ! diagonalize 2p-effective shell model hamiltonian
  CALL zgeev( jobvl, jobvr, np, heff_rhs, lda, ceig_p, vl, ldvl, eigen_vec, ldvr, &
       work1, lwork, rwork, info )
  !   compare spectrum from exact and P-space diagonalization
  WRITE(6,*) 'Compare model space two-body spectrum with exact spectrum:' 
  DO i_p = 1, np
     a1 = REAL( ceig_p(i_p))
     a2 = AIMAG( ceig_p(i_p))
     b1 = REAL( ceig(np+1-i_p) )
     b2 = AIMAG( ceig(np+1-i_p) ) 
     WRITE(6,*) A1, B1
  ENDDO
END SUBROUTINE lee_suzuki2


SUBROUTINE test(np,nq,npq)
  IMPLICIT NONE
  INTEGER,INTENT(IN)::np,nq,npq
  INTEGER::i,j
  COMPLEX*16,DIMENSION(npq,npq)::h
  COMPLEX*16, DIMENSION(npq,npq):: cvec
  COMPLEX*16, DIMENSION(np,np):: cvec_pp
  COMPLEX*16, DIMENSION(nq,np):: cvec_qp
  COMPLEX*16, DIMENSION(npq):: ceig
  COMPLEX*16, DIMENSION(np):: ceig_p
  COMPLEX*16, DIMENSION(np,np)::  heff
  DO i=1,npq
     DO j=i,npq
        h(i,j)=i+j*2
        h(j,i)=i+j*2
     ENDDO
  ENDDO
  CALL vlowkdiag_exact( h, cvec, ceig, npq )
  write(6,*) "test======================="
  DO i=1,npq
     WRITE(6,*) REAL(ceig(i))
  ENDDO
  write(6,*) "======================="
  cvec_pp=0
  cvec_qp=0
  DO i=1,np
     ceig_p(i)=ceig(i)
     DO j=1,np
        cvec_pp(i,j)=cvec(i,j)
     ENDDO
  ENDDO

  DO i=1,nq
     DO j=1,np
        cvec_pp(i,j)=cvec(i+np,j)
     ENDDO
  ENDDO
  call lee_suzuki2(cvec_pp,cvec_qp,ceig_p,np,nq,heff)
END SUBROUTINE test

!
! Lee Suzuki similarity transformation 

SUBROUTINE mylee_suzuki2( cvec_pp, cvec_qp, ceig, np, nq, heff )
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: np, nq
  COMPLEX*16, DIMENSION(np,np), INTENT(IN) :: cvec_pp
  COMPLEX*16, DIMENSION(nq,np), INTENT(IN)  :: cvec_qp
  COMPLEX*16, DIMENSION(np), INTENT(IN) :: ceig
  COMPLEX*16, DIMENSION(np,np) :: temp, omega2, sim1, sim2, omega_2inv, u, u_inv, & 
       cvec_pp_inv
  COMPLEX*16, DIMENSION(nq,np) :: omega
  COMPLEX*16, DIMENSION(np) :: ceig_p, omega2_eig, eigval2
  COMPLEX*16, DIMENSION(np,np) :: eigen_vec, vl, heff_rhs
  COMPLEX*16, DIMENSION(np,np), INTENT(INOUT) ::  heff
  REAL(KIND = 8), DIMENSION(2*np) :: rwork
  COMPLEX*16, DIMENSION(10000) :: work1
  COMPLEX*16 :: d, sum1, temp1(np,np), temp2(np,nq), norm, determinant
  INTEGER :: i_p,j_p, j_q, ii, jj, k, i_q , a_p, a_pp 
  INTEGER :: i, lda, ldb, ldvl, ldvr, info, lwork, ilo , ihi
  CHARACTER*1 :: jobvl, jobvr, balanc, sense
  REAL(KIND = 8), DIMENSION(np) :: scale, rconde, rcondv
  REAL(KIND = 8) :: abnrm
  INTEGER :: ipiv(np) ,j
  REAL(KIND = 8) :: a1, a2, b1, b2

  balanc = 'n';  jobvl = 'n' ;  jobvr = 'v';  sense = 'n';  lda = np
  ldvl = 1;  ldvr = np;  lwork = 10000
  eigen_vec = 0. 
  temp1 = TRANSPOSE(cvec_pp) 
  CALL zgeev( jobvl, jobvr, np, temp1, lda, omega2_eig, vl, ldvl, eigen_vec, ldvr, &
       work1, lwork, rwork, info )
  determinant = PRODUCT(omega2_eig(:))
  !  write(6,*) 'check determinant', determinant

  omega2 = MATMUL(cvec_pp,TRANSPOSE(cvec_pp))
  CALL sqrtmat( omega2, U, U_inv ,np )
  heff=0.
  DO i_p = 1, np
     DO j_p = 1, np
        heff(i_p,j_p) = SUM( cvec_pp(i_p,:)*ceig(:)*cvec_pp(j_p,:) )
     ENDDO
  ENDDO

  heff_rhs =  MATMUL( U_inv, MATMUL( heff,U_inv ) ) 
  ! check if heff is symmetrized:
  heff = 0.
  heff = heff_rhs
  DO i_p = 1, np
     DO j_p = 1, np
        ! make heff manifestly symmetric
        IF ( i_p /= j_p ) THEN
           IF ( ABS(heff(i_p,j_p)- heff(j_p,i_p)) < 1.E-6 ) CYCLE 
           !           WRITE(6,*) 'sym test', heff(i_p,j_p), heff(j_p,i_p)
        ENDIF
     ENDDO
  ENDDO
  ! diagonalize 2p-effective shell model hamiltonian
  CALL zgeev( jobvl, jobvr, np, heff_rhs, lda, ceig_p, vl, ldvl, eigen_vec, ldvr, &
       work1, lwork, rwork, info )
  !   compare spectrum from exact and P-space diagonalization
  WRITE(6,*) 'Compare model space two-body spectrum with exact spectrum:' 
  DO i_p = 1, np
     a1 = REAL( ceig_p(i_p))
     a2 = AIMAG( ceig_p(i_p))
     b1 = REAL( ceig(np+1-i_p) )
     b2 = AIMAG( ceig(np+1-i_p) ) 
     WRITE(6,*) A1, B1
  ENDDO

END SUBROUTINE mylee_suzuki2


!
! eigenvalue sort
! sort cmplx vector real(a(1))<real(a(2)) < ... < real(a(n))
!
SUBROUTINE eigenvalue_sort_cmplx(A, n)
  IMPLICIT NONE
  INTEGER :: i, j, n
  COMPLEX*16, DIMENSION(n), INTENT(INOUT) :: A
  COMPLEX*16 :: temp, temp1
  COMPLEX*16, DIMENSION(n) :: temp2

  DO i = 1, n
     DO j = 1, n
        IF ( ABS( A(i) )  > ABS( A(j) ) ) THEN
           temp = A(i)
           A(i) = A(j) 
           A(j) = temp
        END IF
     END DO
  END DO

END SUBROUTINE eigenvalue_sort_cmplx
!
!fast diagonalizing of complex symmetric matrices
!
SUBROUTINE vlowkdiag_exact( h, cvec, ceig, n )
  IMPLICIT NONE
  INTEGER, INTENT(in) :: n
  COMPLEX*16, DIMENSION(n,n), INTENT(in) :: h
  INTEGER :: k, np
  INTEGER :: p, i,j, i1, kvec, lwork, option, il ,iu, info
  REAL(kind = 8) :: A(n*(n+1)), work(300*n)
  REAL(kind = 8) :: thresh
  COMPLEX*16 :: cvl, cvu
  LOGICAL :: flag(n)  
  COMPLEX*16 :: ceig(n), sum1
  COMPLEX*16 :: cvec(n,n)
  REAL(kind = 8) ::  vl,vu

  lwork = 300*n
  thresh = 30.0
  kvec = n
  option = 4
  info = 1
  i1 = 0
  DO i =  1, n
     DO j =  1, i
        i1 = i1 + 1
        a(i1) = DBLE(h(j,i))
        a(i1+n*(n+1)/2) = AIMAG(h(j,i))
     END DO
  END DO
  CALL cs(n,a,ceig,kvec,cvec,lwork,work,thresh, option,il,iu,cvl,cvu,flag,info)

END SUBROUTINE vlowkdiag_exact
!
!          Obtain the bare potential in oscillator basis
!          Final V is in units of MeV or eV
!          Time consuming part. Parallel version available.
!
SUBROUTINE vlowk_hosc(wave_bra, wave_ket, a,bra,ket,vsum)
  USE wave_functions
  USE constants
  USE relcm_gmatrix
  IMPLICIT NONE
  REAL(KIND = 8), INTENT(INOUT) :: vsum
  INTEGER, INTENT(IN) :: bra, ket
  INTEGER :: i, j
  REAL(KIND = 8), DIMENSION(n_k1,n_k1), INTENT(IN) :: a
  REAL(KIND = 8) :: sum1, sum2, hbarc3
  REAL(KIND = 8), DIMENSION(n_rel), INTENT(IN) :: wave_bra, wave_ket 

  hbarc3=hbarc**3
  sum1=0.
  DO i=1, n_k1
     sum2=0.
     DO j=1, n_k1
        sum2=sum2+wave_ket(j)*a(j,i)
     ENDDO
     sum1=sum1+sum2*wave_bra(i)
  ENDDO
  vsum = sum1*hbarc**3

END SUBROUTINE vlowk_hosc

!
!                 Set Coulomb interaction
!                 in coordinate space, nuclear case, returns in MeV
!
SUBROUTINE coulomb_integral
  USE constants
  USE wave_functions
  USE relcm_gmatrix
  IMPLICIT NONE
  INTEGER :: n1, n2, lr, i, nrel_max
  REAL(DP) :: oscl_r, int_sum, xr, xp, z, factor1, factor2
  REAL(DP), ALLOCATABLE, DIMENSION(:) :: rr, wrr
  REAL(DP) :: cx(0:200)

  oscl_r=oscl*SQRT(2.0D0)        ! Oscillator parameter for relative
  nrel_max = 400
  ALLOCATE ( rr(nrel_max ), wrr(nrel_max ))
  CALL gauss_legendre(0.d0, 20.d0, rr, wrr, nrel_max )
  DO lr = 0, lmax, 1
     DO n1 = 0, nmax, 1
        factor1 = 0.5D0*((n1+lr+2)*LOG(2.D0)+fac(n1)-dfac(2*n1+2*lr+1)-0.5D0*LOG(pi))
        factor1 = EXP(factor1)
        DO n2 = 0, nmax,1
           IF ( n2 == n1) THEN
              factor2 = factor1
           ELSE
              factor2 = 0.5D0*((n2+lr+2)*LOG(2.D0)+fac(n2)-dfac(2*n2+2*lr+1)-0.5D0*LOG(pi))
              factor2 = EXP(factor2)
           ENDIF
           int_sum = 0.D0
           DO i=1,nrel_max
              z= rr(i)/oscl_r
              CALL laguerre_general( n1, lr+0.5D0, z*z, cx )
              xp = cx(n1)*EXP(-z*z*0.5)*(z**lr)
              IF ( n1 == n2) THEN 
                 xr = xp
              ELSE
                 CALL laguerre_general( n2, lr+0.5D0, z*z, cx )
                 xr = cx(n2)*EXP(-z*z*0.5)*(z**lr)
              ENDIF
              int_sum=int_sum+wrr(i)*rr(i)*xp*xr
           ENDDO
           !  Coulomb energy in MeV or eV
           coulomb_relcom(n2, n1, lr) = int_sum*factor1*factor2*1.439965183D0/(oscl_r**3)
        ENDDO
     ENDDO
  ENDDO
  DEALLOCATE ( rr, wrr)

END SUBROUTINE coulomb_integral



!
!           Setup moshinsky brackets for transf
!           rel-cm frame ----> lab frame
!
!
SUBROUTINE mosh_transf(it,ij,gmatrix_configs,lab_to_relcoeff, &
     lab_to_relconf,lab_to_relnumber,max_coeff)
  USE single_particle_orbits
  USE configurations
  IMPLICIT NONE
  TYPE (configuration_descriptor), INTENT(IN) :: gmatrix_configs
  INTEGER, INTENT(IN) :: it, ij, max_coeff 
  INTEGER :: k1,k2,i,j,n1,l1,j1d,n2,l2,j2d,k, number
  REAL(DP), DIMENSION(gmatrix_configs%number_confs,max_coeff), & 
       INTENT(INOUT) :: lab_to_relcoeff
  INTEGER, DIMENSION(gmatrix_configs%number_confs,max_coeff), &
       INTENT(INOUT) :: lab_to_relconf 
  INTEGER, DIMENSION(gmatrix_configs%number_confs), &
       INTENT(INOUT) :: lab_to_relnumber 
  REAL(DP), DIMENSION(max_coeff) :: coeb
  INTEGER, DIMENSION(max_coeff) :: mqnb

  DO k=1,gmatrix_configs%number_confs              
     k2=k*2
     k1=k2-1
     i=gmatrix_configs%config_ab(k1)
     j=gmatrix_configs%config_ab(k2)
     n1=all_orbit%nn(i) 
     l1=all_orbit%ll(i)
     j1d=all_orbit%jj(i)
     n2=all_orbit%nn(j)
     l2=all_orbit%ll(j)
     j2d=all_orbit%jj(j)
     CALL transf_mbs(n1,l1,j1d,n2,l2,j2d,it,ij,mqnb,coeb,max_coeff,number)
     lab_to_relcoeff(k,:)=coeb(:)
     lab_to_relconf(k,:)=mqnb(:)
     lab_to_relnumber(k)=number
  ENDDO

END SUBROUTINE  mosh_transf

SUBROUTINE mosh_transf_t(it,itz,ij,gmatrix_configs,lab_to_relcoeff, &
     lab_to_relconf,lab_to_relnumber,max_coeff)
  USE single_particle_orbits
  USE configurations
  IMPLICIT NONE
  TYPE (configuration_descriptor), INTENT(IN) :: gmatrix_configs
  INTEGER, INTENT(IN) :: it, itz, ij, max_coeff 
  INTEGER :: k1,k2,i,j,n1,l1,j1d,n2,l2,j2d,k, number
  REAL(DP) ::fac,fac1,delta
  REAL(DP), DIMENSION(gmatrix_configs%number_confs,max_coeff), & 
       INTENT(INOUT) :: lab_to_relcoeff
  INTEGER, DIMENSION(gmatrix_configs%number_confs,max_coeff), &
       INTENT(INOUT) :: lab_to_relconf 
  INTEGER, DIMENSION(gmatrix_configs%number_confs), &
       INTENT(INOUT) :: lab_to_relnumber 
  REAL(DP), DIMENSION(max_coeff) :: coeb
  INTEGER, DIMENSION(max_coeff) :: mqnb

  DO k=1,gmatrix_configs%number_confs              
     k2=k*2
     k1=k2-1
     i=gmatrix_configs%config_ab(k1)
     j=gmatrix_configs%config_ab(k2)
     n1=all_orbit%nn(i) 
     l1=all_orbit%ll(i)
     j1d=all_orbit%jj(i)
     n2=all_orbit%nn(j)
     l2=all_orbit%ll(j)
     j2d=all_orbit%jj(j)
     CALL transf_mbs_t(n1,l1,j1d,n2,l2,j2d,it,itz,ij,mqnb,coeb,max_coeff,number)
     lab_to_relcoeff(k,:)=coeb(:)
     lab_to_relconf(k,:)=mqnb(:)
     lab_to_relnumber(k)=number
!     IF (itz /=0) THEN
!        CALL transf_mbs(n1,l1,j1d,n2,l2,j2d,itz,ij,mqnb,coeb1,max_coeff,number)
!        lab_to_relcoeff(k,:)=coeb1(:)
!        lab_to_relconf(k,:)=mqnb(:)
!        lab_to_relnumber(k)=number
!     ELSE
!        CALL transf_mbs(n1,l1,j1d,n2,l2,j2d,itz,ij,mqnb,coeb1,max_coeff,number)
!        CALL transf_mbs(n2,l2,j2d,n1,l1,j1d,itz,ij,mqnb,coeb2,max_coeff,number)
!        fac = 1.D0
!        fac1=1.D0
!        IF (MOD(ij+(j1d+j2d)/2+1,2) /= 0 ) fac1= -1.D0
!        !    if 1==2,ij+it must be odd
!        IF (it == 0) THEN
!           fac=1.D0/(SQRT(2.D0)*SQRT(1.d0+delta(n1,n2)*delta(l1,l2)*delta(j1d,j2d)))
!           lab_to_relcoeff(k,:)=fac*(coeb1(:)-fac1*coeb2(:))
!           lab_to_relconf(k,:)=mqnb(:)
!           lab_to_relnumber(k)=number
!
!        ELSE
!           fac=1.D0/(SQRT(2.D0)*SQRT(1.d0+delta(n1,n2)*delta(l1,l2)*delta(j1d,j2d)))
!           lab_to_relcoeff(k,:)=fac*(coeb1(:)+fac1*coeb2(:))
!           lab_to_relconf(k,:)=mqnb(:)
!           lab_to_relnumber(k)=number
!        END IF
!     END IF
  ENDDO

END SUBROUTINE  mosh_transf_t



SUBROUTINE transf_mbs_t(n1,l1,j1d,n2,l2,j2d,it,itz,ij,mqnb,coeb,max_coeff,number)
  USE single_particle_orbits
  USE partial_waves
  USE ang_mom_functions
  USE configurations
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: n1,l1,j1d,n2,l2,j2d,it,itz,ij,max_coeff
  INTEGER, DIMENSION(max_coeff), INTENT(OUT) :: mqnb
  REAL(DP), DIMENSION(max_coeff), INTENT(OUT) :: coeb
  INTEGER :: ld, l, n, lcd, lambd, isign, lambdd,  a, b, k1, k2, &
       l2d, l1d, nemx, ijd, l_min, l_max, &
       jlsd, nch, isd, jls, ispin,nc, lc, nconfs
  INTEGER, INTENT(OUT) :: number
  REAL(DP) :: bk, w6, w9, s , co, cst, delta
  LOGICAL triag

  mqnb=0 ; coeb=0.
  l1d=l1*2 ;   l2d=l2*2 ;   ijd=ij*2 ;   nemx=2*n1+l1+2*n2+l2
  number=0
  !     loop over possible partial wave channels which specify J, l, S and Tz
  DO nch=1,no_channels
     IF ( itz /= iso(nch) ) CYCLE
     l_min=orb_lrel_min(nch)
     l_max=orb_lrel_max(nch)
     ispin=spin_rel(nch) ; isd=2*ispin
     jls=jang_rel(nch) ; jlsd=2*jls
     !     loop over rel and cm configs specifying n, l, N and L
     DO nconfs =1, relcm_conf%nconfs_relcm(nch)
        k2=nconfs*2
        k1=k2-1
        b=relcm_conf%relcm_ab(nch,k2)
        a=relcm_conf%relcm_ab(nch,k1)
        n=relcm_sp_data%nrel(a)
        l=relcm_sp_data%lrel(a)
        IF ( (l /= l_max) .AND. ( l /= l_min ) ) CYCLE
        nc=relcm_sp_data%nrel(b)
        lc=relcm_sp_data%lrel(b)
        !     pauli test for identical particles in partial waves
        IF((-1)**(l+ispin+it) > 0) CYCLE
        cst = SQRT(2.D0)/SQRT( 1.d0+ delta(n1,n2)*delta(l1,l2)*delta(j1d,j2d) )
        IF(triag(jls,ij,lc)) CYCLE
        IF((-1)**(l1+l2+l+lc) < 0) CYCLE
        IF((2*n+l+2*nc+lc) /= nemx) CYCLE
        IF((2*n+l+2*nc+lc) > nlmax) CYCLE
        IF(triag(jls,l,ispin)) CYCLE
        ld=l*2 ; lcd=lc*2
        co=0.D0
        DO lambd = ABS(l1-l2), l1+l2
           lambdd=lambd*2
           IF(triag(lambd,l,lc)) CYCLE
           IF(triag(ij,lambd,ispin)) CYCLE
           bk= gmosh(n,l,nc,lc,n1,l1,n2,l2,lambd,1.D0)
           !           IF( ABS (bk) == 0.D0) CYCLE
           w9=snj(l1d,1,j1d,l2d,1,j2d,lambdd,isd,ijd)
           !           IF ( ABS( w9) == 0.D0) CYCLE
           w9=w9*SQRT((j1d+1.)*(j2d+1.)*(lambdd+1.)*(isd+1.))
           isign=(lcd+ld+ijd+isd)/2
           s=1.-2.*MOD(isign,2)
           w6=s*sjs(lcd,ld,lambdd,isd,ijd,jlsd)
           w6=w6*SQRT(FLOAT((lambdd+1)*(jlsd+1)))
           co=co+cst*bk*w6*w9*(-1.)**(l+lambd+jls+ij)
        ENDDO
        IF(DABS(co) <= 1.D-8) co = 0.0_dp!THEN
        number=number+1
        mqnb(number)=nconfs*1000+nch
        coeb(number)=co
        !        ENDIF
     ENDDO
  ENDDO

END SUBROUTINE transf_mbs_t
!
!           Explicit evaluation of transf coeff, refer to prc,19,3(1979)
!
SUBROUTINE transf_mbs(n1,l1,j1d,n2,l2,j2d,itz,ij,mqnb,coeb,max_coeff,number)
  USE single_particle_orbits
  USE partial_waves
  USE ang_mom_functions
  USE configurations
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: n1,l1,j1d,n2,l2,j2d,itz,ij,max_coeff
  INTEGER, DIMENSION(max_coeff), INTENT(OUT) :: mqnb
  REAL(DP), DIMENSION(max_coeff), INTENT(OUT) :: coeb
  INTEGER :: ld, l, n, lcd, lambd, isign, lambdd,  a, b, k1, k2, &
       l2d, l1d, nemx, ijd, l_min, l_max, &
       jlsd, nch, isd, jls, ispin,nc, lc, nconfs
  INTEGER, INTENT(OUT) :: number
  REAL(DP) :: bk, w6, w9, s , co, cst, delta
  LOGICAL triag

  mqnb=0 ; coeb=0.
  l1d=l1*2 ;   l2d=l2*2 ;   ijd=ij*2 ;   nemx=2*n1+l1+2*n2+l2
  number=0
  !     loop over possible partial wave channels which specify J, l, S and Tz
  DO nch=1,no_channels
     IF ( itz /= iso(nch) ) CYCLE
     l_min=orb_lrel_min(nch)
     l_max=orb_lrel_max(nch)
     ispin=spin_rel(nch) ; isd=2*ispin
     jls=jang_rel(nch) ; jlsd=2*jls
     !     loop over rel and cm configs specifying n, l, N and L
     DO nconfs =1, relcm_conf%nconfs_relcm(nch)
        k2=nconfs*2
        k1=k2-1
        b=relcm_conf%relcm_ab(nch,k2)
        a=relcm_conf%relcm_ab(nch,k1)
        n=relcm_sp_data%nrel(a)
        l=relcm_sp_data%lrel(a)
        IF ( (l /= l_max) .AND. ( l /= l_min ) ) CYCLE
        nc=relcm_sp_data%nrel(b)
        lc=relcm_sp_data%lrel(b)
        !     pauli test for identical particles in partial waves
        IF((ABS(itz) == 1).AND. ((-1)**(l+ispin+1) > 0)) CYCLE
        !     get right sign from isospin clebsch-gordan coeff for the partial waves
        !     since we have for Tz=0 always a coupling order | (pn) J >
        !     this means that the T=0 contrib goes like  -(1-(-)^(l+spin))/sqrt(2)
        !     and l+s have to be odd       
        cst = 1.0_dp
        IF ( itz /= 0 ) THEN
           cst = SQRT(2.D0)/SQRT( 1.d0+ delta(n1,n2)*delta(l1,l2)*delta(j1d,j2d) )
        ELSEIF ( (itz == 0) ) THEN
           IF ( (MOD(l+ispin,2) /= 0 ) )  cst= -1.0_dp
        ENDIF
        IF(triag(jls,ij,lc)) CYCLE
        IF((-1)**(l1+l2+l+lc) < 0) CYCLE
        IF((2*n+l+2*nc+lc) /= nemx) CYCLE
        IF((2*n+l+2*nc+lc) > nlmax) CYCLE
        IF(triag(jls,l,ispin)) CYCLE
        ld=l*2 ; lcd=lc*2
        co=0.D0
        DO lambd = ABS(l1-l2), l1+l2
           lambdd=lambd*2
           IF(triag(lambd,l,lc)) CYCLE
           IF(triag(ij,lambd,ispin)) CYCLE
           bk= gmosh(n,l,nc,lc,n1,l1,n2,l2,lambd,1.D0)
!           IF( ABS (bk) == 0.D0) CYCLE
           w9=snj(l1d,1,j1d,l2d,1,j2d,lambdd,isd,ijd)
!           IF ( ABS( w9) == 0.D0) CYCLE
           w9=w9*SQRT((j1d+1.)*(j2d+1.)*(lambdd+1.)*(isd+1.))
           isign=(lcd+ld+ijd+isd)/2
           s=1.-2.*MOD(isign,2)
           w6=s*sjs(lcd,ld,lambdd,isd,ijd,jlsd)
           w6=w6*SQRT(FLOAT((lambdd+1)*(jlsd+1)))
           co=co+cst*bk*w6*w9*(-1.)**(l+lambd+jls+ij)
        ENDDO
        IF(DABS(co) <= 1.D-8) co = 0.0_dp!THEN
           number=number+1
           mqnb(number)=nconfs*1000+nch
           coeb(number)=co
!        ENDIF
     ENDDO
  ENDDO

END SUBROUTINE transf_mbs
!
!           Finds total number of transformation  coeffs 
!           for given J, parity and isospin projection
!
SUBROUTINE find_ncoeffs(it,ij,this,max_coeff)
  USE single_particle_orbits
  USE configurations
  IMPLICIT NONE
  TYPE (configuration_descriptor), INTENT(IN) :: this
  INTEGER, INTENT(IN) :: it, ij 
  INTEGER :: k1,k2,i,j,n1,l1,n2,l2,number, k
  INTEGER, INTENT(OUT) :: max_coeff
  max_coeff=0
  DO k=1, this%number_confs              
     k2=k*2
     k1=k2-1
     i=this%config_ab(k1)
     j=this%config_ab(k2)
     n1=all_orbit%nn(i) 
     l1=all_orbit%ll(i)
     n2=all_orbit%nn(j)
     l2=all_orbit%ll(j)
     CALL count_trans_mbs(n1,l1,n2,l2,it,ij,number)
     IF ( max_coeff < number ) max_coeff = number
  ENDDO

END SUBROUTINE find_ncoeffs
!
!           counts max possible number  of transf coeff
!
SUBROUTINE count_trans_mbs(n1,l1,n2,l2,itz,ij,number)
  USE single_particle_orbits
  USE partial_waves
  USE configurations
  IMPLICIT NONE
  INTEGER, INTENT(OUT) :: number
  INTEGER, INTENT(IN) :: n1,l1,n2,l2,itz,ij
  INTEGER :: l, n, nemx, l_min, l_max, nconfs , nch, jls, & 
       ispin, nc, lc, a, b,  k1, k2 
  LOGICAL triag

  nemx=2*n1+l1+2*n2+l2
  number=0
  DO nch=1,no_channels
     IF ( itz /= iso(nch) ) CYCLE
     l_min=orb_lrel_min(nch)
     l_max=orb_lrel_max(nch)
     ispin=spin_rel(nch)
     jls=jang_rel(nch)
     DO nconfs =1, relcm_conf%nconfs_relcm(nch)
        k2=nconfs*2
        k1=k2-1
        b=relcm_conf%relcm_ab(nch,k2)
        a=relcm_conf%relcm_ab(nch,k1)
        n=relcm_sp_data%nrel(a)
        l=relcm_sp_data%lrel(a)
        IF ( (l /= l_max) .AND. ( l /= l_min ) ) CYCLE
        nc=relcm_sp_data%nrel(b)
        lc=relcm_sp_data%lrel(b)
        IF((ABS(itz)==1).AND.((-1)**(l+ispin+1) > 0)) CYCLE
        IF(triag(jls,ij,lc)) CYCLE
        IF((-1)**(l1+l2+l+lc) < 0) CYCLE
        IF((2*n+l+2*nc+lc) /= nemx) CYCLE
        IF((2*n+l+2*nc+lc) > nlmax) CYCLE
        IF(triag(jls,l,ispin)) CYCLE
        number=number+1
     ENDDO
  ENDDO

END SUBROUTINE count_trans_mbs

!
!                 Set up h.o. wf for rel system
!
SUBROUTINE nocoreho_wfunction
  USE constants
  USE wave_functions
  IMPLICIT NONE
  INTEGER :: n, l, i, j
  REAL(KIND = 8)  :: cx(0:200), factor, z_rel, xp, ph, oscl_r, sum_rel, contrib

  oscl_r=oscl*SQRT(2.)        ! Oscillator parameter for relative
  DO n=0,nmax
     ph=(-1.D0)**n
     DO l=0,lmax
        sum_rel=0.
        factor = 0.5D0*((n+l+2)*LOG(2.D0)+fac(n)-dfac(2*n+2*l+1)-0.5D0*LOG(pi))
        factor = EXP(factor)
        DO i=1,n_rel
           z_rel= ra(i)*ra(i)*oscl_r*oscl_r
           CALL laguerre_general( n, l+0.5D0, z_rel, cx )
           xp = EXP(-z_rel*0.5)*((ra(i)*oscl_r)**l)*cx(n)
           rnlr(i,l,n) = xp*(wra(i)*(ra(i)**2))*ph*factor*(oscl_r**(1.5D0))     ! rel wf
           contrib = xp*ph*factor*(oscl_r**(1.5D0)) 
           sum_rel=sum_rel+ wra(i)*(contrib*ra(i))**2
        ENDDO
        WRITE(6,*) 'Norm rel ho wf n,l : ', n, l, sum_rel
     ENDDO
  ENDDO

END SUBROUTINE nocoreho_wfunction

SUBROUTINE nocorepotential_interface(ncoup,vkk,ichan)
  USE wave_functions
  USE relcm_gmatrix
  USE partial_waves
  USE constants
  IMPLICIT NONE
  INTEGER :: i,j, ncoup, k, l, jxx, ichan, inn, isospin_tz, spin, &
       n1, ix, iy
  REAL(KIND=8)  :: v,xmev,ymev, c, q
  CHARACTER (LEN=4) :: label
  COMMON /cnn/ inn
  COMMON /cpot/ v(6),xmev,ymev
  common /cpts/   q(97),c,n1,ix,iy
  COMMON /cstate/ jxx, heform, sing, trip, coup, endep, label
  LOGICAL :: sing, trip, coup, heform, endep
  REAL(KIND=8), DIMENSION(ncoup*n_rel,ncoup*n_rel), INTENT(OUT) :: vkk
  REAL(KIND=8) :: v00, v11, v12, v21, v22, fmult, besl, bs, b11, b12, b21, b22
  REAL(KIND=8) :: xb(500), wb(500), v8(500,2,2), temp(500,2,2)
  INTEGER :: j1,l1,l2,s1,isot,m,nt,tz1,tz2,lpot

  heform=.FALSE.
  jxx=jang_rel(ichan)
  sing=spin_rel(ichan) == 0
  trip=spin_rel(ichan) == 1
  coup=(ncoup == 2) 
  isospin_tz=iso(ichan)
  spin = spin_rel(ichan)
  SELECT CASE (type_of_pot)
  CASE('Idaho-A')
     inn = 1
  CASE('Idaho-B')
     inn = 2
  CASE('n3lo')
     SELECT CASE ( isospin_tz)
     CASE (-1)
        IF ( (csb == 'csb').AND.(cib == 'cib'))       inn = 1    !  pp case = 1
        IF ( (csb == 'no-csb').AND.(cib == 'cib'))    inn = 3    !  pp case = nn case
        IF ( (cib == 'no-cib').AND.(csb == 'no-csb')) inn = 2    !  pp case = pn case
        tz1 = 1; tz2 = 1
     CASE (0)
        inn = 2    !  pn case = 2  if all inn=2, no CSB or ISB
        tz1 = 1; tz2 = -1
     CASE ( 1)
        IF ( (csb == 'csb').AND.(cib == 'cib'))       inn = 3    !  nn case = 3
        IF ( (csb == 'no-csb').AND.(cib == 'cib'))    inn = 3    !  nn case = pp case 
        IF ( (cib == 'no-cib').AND.(csb == 'no-csb')) inn = 2    !  nn case = pn case
        tz1 = -1; tz2 = -1
     END SELECT
!!$  CASE('n3lo3b')
!!$     SELECT CASE ( isospin_tz)
!!$     CASE (-1)
!!$        IF ( (csb == 'csb').AND.(cib == 'cib'))       inn = 1    !  pp case = 1
!!$        IF ( (csb == 'no-csb').AND.(cib == 'cib'))    inn = 3    !  pp case = nn case
!!$        IF ( (cib == 'no-cib').AND.(csb == 'no-csb')) inn = 2    !  pp case = pn case
!!$        tz1 = 1; tz2 = 1
!!$     CASE (0)
!!$        inn = 2    !  pn case = 2  if all inn=2, no CSB or ISB
!!$        tz1 = 1; tz2 = -1
!!$     CASE ( 1)
!!$        IF ( (csb == 'csb').AND.(cib == 'cib'))       inn = 3    !  nn case = 3
!!$        IF ( (csb == 'no-csb').AND.(cib == 'cib'))    inn = 3    !  nn case = pp case 
!!$        IF ( (cib == 'no-cib').AND.(csb == 'no-csb')) inn = 2    !  nn case = pn case
!!$        tz1 = -1; tz2 = -1
!!$     END SELECT
  CASE('CD-bonn')
     SELECT CASE ( isospin_tz)
     CASE (-1)
        IF ( (csb == 'csb').AND.(cib == 'cib'))       inn = 1    
        IF ( (csb == 'no-csb').AND.(cib == 'cib'))    inn = 3    
        IF ( (cib == 'no-cib').AND.(csb == 'no-csb')) inn = 2    
        tz1 = 1; tz2 = 1
     CASE (0)
        inn = 2    !  pn case = 2  if all inn=2, no CSB or ISB
        tz1 = 1; tz2 = -1
     CASE ( 1)
        IF ( (csb == 'csb').AND.(cib == 'cib'))       inn = 3    
        IF ( (csb == 'no-csb').AND.(cib == 'cib'))    inn = 3    
        IF ( (cib == 'no-cib').AND.(csb == 'no-csb')) inn = 2    
        tz1 = -1; tz2 = -1
     END SELECT
  END SELECT
  IF ( (type_of_pot == 'argonnev8').OR.(type_of_pot == 'argonnev18').OR.   &
       (type_of_pot == 'OPEP').OR.(type_of_pot == 'Tensorinteraction')   &
       .OR.(type_of_pot == 'LSinteraction').OR.(type_of_pot == 'MalflietTjon') ) THEN
     SELECT CASE ( isospin_tz)
     CASE (-1)
        fmult = 2.D0/ACOS(-1.D0)/p_mass(-1)/hbarc
        tz1 = 1; tz2 = 1
     CASE (0)
        tz1 = 1; tz2 = -1
        fmult = 2.D0/ACOS(-1.D0)/p_mass(0)/hbarc
     CASE ( 1)
        tz1 = -1; tz2 = -1
        fmult = 2.D0/ACOS(-1.D0)/p_mass(1)/hbarc
     END SELECT
     lpot = 1  ! switch for argonne potentials, lpot 1 means v18
     IF (type_of_pot == 'argonnev8') lpot = 2   ! lpot = 2 means v8
     IF (type_of_pot == 'OPEP') lpot = 3   ! lpot = 3 means OPEP from V18
     IF (type_of_pot == 'Tensorinteraction') lpot = 4   ! lpot = 4 means full tensor from V18
     IF (type_of_pot == 'LSinteraction') lpot = 5   ! lpot = 5 means full LS from V18
     IF (type_of_pot == 'MalflietTjon') lpot = 6   ! lpot = 6 Malfliet-Tjon
     !  sett mesh points in r-space
     nt = 300
     CALL gauss_legendre(0.D0,10.D0,xb,wb,nt)
     v8 = 0.
     j1 = jang_rel(ichan); s1 = spin_rel(ichan)
     IF ( ncoup == 1) THEN
        l1 = jang_rel(ichan)
        l2 = jang_rel(ichan)
     ELSEIF ( ncoup == 2) THEN
        l1 = jang_rel(ichan)-1
        l2 = jang_rel(ichan)+1
     ENDIF
     IF ( ABS(isospin_tz) == 1) THEN
        isot = 1
        CALL argonne_pots(j1,l1,s1,isot,tz1,tz2,v8,xb,nt,lpot)
     ELSEIF ( isospin_tz == 0) THEN
        v8 = 0.
        !  now check partial waves in good isospin formalism
        DO isot = 0, 1
           IF(MOD(l1+s1+isot,2) == 0) CYCLE
           temp = 0.
           CALL argonne_pots(j1,l1,s1,isot,tz1,tz2,temp,xb,nt,lpot)
           v8 = v8+temp
        ENDDO
     ENDIF
     DO i=1,n_rel
        IF (ncoup == 2)   k=i+n_rel 
        DO j=1,n_rel
           IF (ncoup == 2) l=j+n_rel 
           IF (ncoup == 1 ) THEN
              v00 = 0.
              DO m=1,nt
                 bs=besl(xb(m)*ra(i),l1)*besl(xb(m)*ra(j),l1)
                 v00=v00+bs*v8(m,1,1)*wb(m)*xb(m)**2
              ENDDO
              vkk(j,i) = v00*fmult
           ELSEIF (ncoup ==2  ) THEN 
              v11=0.d0
              v12=0.d0
              v21=0.d0
              v22=0.d0
              DO m=1,nt
                 c=wb(m)*xb(m)**2
                 b11=besl(xb(m)*ra(i),l1)
                 b12=besl(xb(m)*ra(i),l2)
                 b21=besl(xb(m)*ra(j),l1)
                 b22=besl(xb(m)*ra(j),l2)
                 v11=v11+c*b11*b21*v8(m,1,1)
                 v12=v12+c*b11*b22*v8(m,1,2)
                 v21=v21+c*b12*b21*v8(m,2,1)
                 v22=v22+c*b12*b22*v8(m,2,2)
              ENDDO
              vkk(j,i)=v11*fmult
              vkk(j,k)=v21*fmult
              vkk(l,i)=v12*fmult
              vkk(l,k)= v22*fmult
           ENDIF
        ENDDO
     ENDDO
  ELSE
     DO i=1,n_rel
        xmev=ra(i)*hbarc
        IF (ncoup == 2)   k=i+n_rel 
        DO j=1,n_rel
           ymev=ra(j)*hbarc
           IF (ncoup == 2) l=j+n_rel 
           SELECT CASE (type_of_pot)
           CASE('CD-bonn')
              CALL cdbonn
           CASE('Idaho-A')
              CALL idaho
           CASE('Idaho-B')
              CALL idaho
           CASE('n3lo')
              CALL n3lo
!!$           CASE('n3lo3b')
!!$              CALL n3lo3b
           END SELECT
           IF (sing ) THEN
              vkk(j,i)=v(1)
           ELSEIF ((trip).AND.(.NOT.coup )  ) THEN 
              vkk(j,i) = v(2)
           ELSEIF (coup ) THEN
              vkk(j,i)=v(4)
              vkk(j,k)=v(5)
              vkk(l,i)=v(6)
              vkk(l,k)= v(3)
           ENDIF
        ENDDO
     ENDDO
  ENDIF

END SUBROUTINE nocorepotential_interface
